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We present a hydrodynamic lattice gas model for two-dimensional flows on curved surfaces with 
dynamical geometry. This model is an extension to two dimensions of the dynamical geometry 
lattice gas model previously studied in one-dimension pj[3]. We expand upon a variation of the two- 
dimensional flat space FHP model created by Frisch, Hasslacher and Pomeau, and independently 
by Wolfram [HE], an d modified by Boghosian, Love, and Meyer in [6]. We define a hydrodynamic 
lattice gas model on an arbitrary triangulation, whose flat space limit is the FHP model. Rules that 
change the geometry are constructed using the Pachner moves, which alter the triangulation but 
not the topology [7j. We present results on the growth of the number of triangles as a function of 
time. Simulations show that the number of triangles lattice grows with time as 1 3 , in agreement a 
mean field prediction. We also present preliminary results on the distribution of curvature over a 
typical triangulation for these simulations. 



I. INTRODUCTION 

Lattice-gas automata (LGA) models for fluids date 
from the sixties, when Kadanoff and Swift and Hardy, 
de Pazzis and Pomeau introduced the first such mod- 
els [HI [H]- Both of these models use a two-dimensional 
Cartesian lattice and are anisotropic. Since simple fluids 
are isotropic, these models are not capable of reproducing 
hydrodynamics. This problem was solved in 1986 when 
Frish, Hasslacher, and Pomeau, and independently Wol- 
fram, introduced an isotropic model (the FHP model) 
using a triangular lattice. They demonstrated that an 
LGA models the Navier-Stokes equations in flat two- 
dimensional space [H [S] . 

All LGA methods are characterized by phases of prop- 
agation and collision of particles that move on a lattice. 
During the propagation phase, particles move from site 
to site on the lattice, while during the collision phase 
the particles rearrange themselves amongst the vectors 
at each site (see Figure [1} . Before we discuss the FHP 
rules in detail, it is important to note that the rules that 
govern these models are not meant to replicate the phys- 
ical world on a small scale; the Navier-Stokes equations 
emerge from the FHP rules on the macroscopic scale 
for large lattice sizes and spatial or ensemble averaging. 
The microscopic rules are only required to conserve to- 
tal momentum, particle number, and energy. Addition- 
ally, the lattice must be sufficiently symmetric to yield 
an isotropic pressure tensor. 

Many 2-D situations of physical interest use a Eu- 
clidean plane as the underlying geometry, hence "lattice" 
gases, in which the model is constructed on a translation 
invariant discretization of Euclidean space. However sit- 
uations exist, such as atmospheric flow, the experiments 
of Seychelles [TU] , or surface flows in interfaces embedded 
in fluid mixtures, in which a discretization of a sphere 
or other surface in which the geometry is non-Euclidean 
may be more appropriate. In such geometries, the angles 
of a triangle need not sum to tt. We may specialize to 



simplicial complexes made up of equilateral triangles, as 
any 2-D surface may be discretized in this way [7]. In 
this case the geometry is defined locally by the number 
of triangles meeting at each grid point. If six triangles 
meet, the geometry is locally flat. If fewer than six trian- 
gles meet, the geometry has positive local curvature. If 
more than six triangles meet, the geometry has negative 
local curvature. If the properties of the triangulation, 
including the local curvature, are allowed to change we 
call the geometry dynamical. 

There are many situations in physics in which geom- 
etry takes on a dynamical role. Perhaps the most fun- 
damental is in Einstein's general theory of relativity, in 
which the idea of motion along geodesies in a Rieman- 
nian manifold supervenes Newtonian ideas of accelera- 
tion due to forces [TT] . In the Regge treatment of general 
relativity |12] and the causal dynamical triangulations 
approach to quantum gravity [13l [14] these Riemannian 
manifolds are replaced by simplicial complexes. The sta- 
tistical mechanics and growth dynamics of random sur- 
faces has been much studied for both 1-D interfaces [TH] 
and 2-D surfaces [16TH9j . In spite of their origin in very 
different physical systems, the common language of dis- 
cretized surfaces can be informative. For example, the 
crumpling transition of membranes [19] also occurs in 
Euclidean approaches to simplical quantum gravity [13j . 

In this paper we present a hydrodynamic lattice gas 
model for two-dimensional flows on curved surfaces with 
dynamical geometry. We extend a variation of the FHP 
model to arbitrary equilateral triangulations. We allow 
the geometry so defined to become dynamical by apply- 
ing the Pachner moves contingent on the particle con- 
tent. The restriction of time-reversibility is used to re- 
strict the rule space, as in the one dimensional version 
of this model [TJ-[3]. We present a mean-field prediction 
and simulation results for the growth of the lattice as a 
function of time, and give preliminary results on the dis- 
tribution of curvature on the triangulations generated by 
these simulations. We close the paper with some conclu- 
sions and directions for future work. 
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II. THE FHP MODEL 

In the FHP lattice gas automata the particles move 
on a triangular lattice. At each lattice site there are six 
lattice vectors. Each vector can be occupied by at most 
one particle - the model has an exclusion principle. The 
vector which a particle occupies defines its velocity. In 
Figure [T] for example, the site is occupied by a single 
particle moving to the right. The state of a particular 
site is given by its particle content. Each vector at each 
site can have two states: occupied or unoccupied. This 
gives a total of 2 6 = 64 states per site. 

To facilitate the generalization from a Euclidean lattice 
to an arbitrary triangulation, we would like to regard our 
sites as triangles rather than a single point. We there- 
fore enclose the site in a triangle and move the vectors to 
the edges of the triangles, as demonstrated by Figure [T] 
This "inflated" site is equivalent to the point site used by 
Frisch, Hasslacher and Pomeau. This modification of the 
FHP model was originally proposed by Boghosian, Love, 
and Meyer |6j . Those authors proposed collisions on the 
edges of the triangles, where four vectors meet and in- 
troduced the possibility of having a rest particle of mass 
two in the model for a total of five bits per state. This 
model was analyzed by a grouping of triplets of trian- 
gle edges sites followed by a Chapman-Enskog expansion 
which yielded isotropic fluid equations. 

However, running a channel flow simulation using their 
proposed model produced the image in Figure [4] Note 
that the structure of the lattice is evident on a macro- 
scopic scale in the figure. This is due to a spurious 
conserved momentum in collisions at the vertices of the 
Kagome lattice. No momentum is transferred between 
separate lines of the lattice, so momentum is conserved 
in three directions in two-dimensional space. This leads 
to unphysical flows, an example of which is shown in Fig- 
ure [U We therefore redefine the FHP model with colli- 
sions occurring on the faces, rather than the edges, of the 
triangle. It should be noted that the Chapman-Enskog 
analysis presented in [6] remains valid for a model, such 
as the one we present here, in which momentum is ex- 
changed by collisions among all lattice directions. 

The rules used for fixed geometry in the variant of the 
FHP model we study are shown graphically in Figures [2] 
and [3] If two particles enter a site with opposite veloci- 
ties, as in Figure|2j they flip to either of the other lines of 
the lattice with equal probability. If three particles enter 
a site such that their total momentum sums to zero, in 
other words, there is a particle occupying every other vec- 
tor, the particles switch from the occupied vectors to the 
unoccupied vectors, shown in Figure [3j This three-body 
collision breaks the separate conservation of momentum 
along each line of the lattice. This is required because 
extra conserved quantities lead to incorrect macroscopic 
behavior. If particles enter in any other configuration, 
they are simply allowed to propagate as usual to the next 
site along their geodesic [I]. 




FIG. 1: An FHP lattice site has 6 possible velocities labeled 
through 5, each of which represents the velocity of a particle. 
Each vector can hold at most one particle, so that each site 
has 2 6 = 64 states. The traditional representation of a site 
in the FHP model is the star shown on the left. By moving 
the vectors to the edges of a triangle, as shown in the center 
picture, we convert the site from a single point to the face 
of a triangle. These two sites are equivalent. If we remove 
the arrow heads from the vectors, we produce the notation- 
ally convenient right hand figure. We refer to the conversion 
between the star and the triangle as inflation. 




FIG. 2: A two-particle collision at a site. The particles switch 
with equal probability to one of the other two directions of 
the lattice. This rule applies to any two particles entering a 
site with opposite velocities. 



III. LATTICE GASES ON CURVED SURFACES 

We now generalize the FHP model to arbitrary equilat- 
eral triangulations. It is known that any manifold can be 
approximated arbitrarily closely by a tiling of equilateral 
triangles [7]. This allows us to triangulate any surface 
and regard each face as an inflated FHP site. In the spe- 
cial case of flat space the vectors in the array of inflated 
sites (Figure [I]) create a tiling of Stars of David, along 
the lines of which the particles can move. This lattice is 
known as the Kagome lattice. Figure [5] shows a triangu- 
lation of a cylinder and an icosahedron where the nodes 
of the triangulation are shown in white and the nodes of 
the Kagome lattice are in red. These images were gener- 
ated with visual python [2D] . We now describe the rules 
which couple the particles to the triangulation and allow 
the geometry to become dynamical. 

To allow the geometry to become dynamical, we em- 
ploy the Pachner moves. A sequence of Pachner moves 
cannot change the topology of a manifold, but it can take 
the manifold from one triangulation to another: a torus 
can morph into another toroidal geometry such as a cof- 
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FIG. 3: A three-particle collision at a site. The particles 
switch to the unoccupied vectors of the site. This move breaks 
the separate conservation of momentum along each direction 
of the lattice. 




FIG. 4: Channel flow simulation with a barrier using the 
model defined in [£]. Arrows indicate the velocity field, and 
the color scale indicates vorticity. In this formulation of the 
model collisions occur on the edges, rather than the faces, of 
the triangulation. This leads to separate conservation of mo- 
mentum along each lattice direction and hence an anisotropic 
model. This anisotropy is evident in this simulation as the 
structure of the lattice is visible in the flow. 




FIG. 5: A Kagome lattice on a cylinder (left) and on an icosa- 
hedron (right). The nodes of the triangulation are depicted as 
white spheres. Any triangulation of a two-dimensional mani- 
fold can be tiled with a Kagome lattice. Particles are allowed 
to move along the lines of the Kagome lattice, and are shown 
in black on the left. These images were generated with visual 
python [20] . 




FIG. 6: Pachner Moves. A sequence of Pachner moves can 
connect any pair of triangulations of a manifold, but cannot 
change the topology. The two-to-two move (left) changes the 
orientation of two triangles as show above, and the one-to- 
three move (right) replaces one triangle with three, creating 
a tetrahedron, or vice versa. 



fee mug, but it can not morph into a sphere [TJ. The state 
of the system with static and flat geometry is specified 
by the particle content of the sites alone, the state of the 
system with dynamical geometry specified by both the 
particle content and the geometry of the triangulation. 

There are two Pachner moves for two dimensional tri- 
angulations: a two-to-two move, where the number of tri- 
angles is unchanged, and a one-to-three or three-to-one 
move that increases or decreases the number of triangles 
by two. We call the three-to-one and the one-to-three 
move addition and deletion, respectively, because they 
add or subtract a tetrahedron from the surface. 

In three dimensions the two-to-two Pachner move (Fig- 
ure^ left) is not isometrically embeddable in general. If 
two triangles are removed, turned, and replaced in the 
triangulation, they will not fit unless the dihedral an- 
gle between the original pair of triangles was that of the 
tetrahedron. This is unsurprising since two-dimensional 
manifolds are not generically embeddable in three dimen- 
sions [21 . However, the one-to-three move, (Figure |6j 
right), is generically immersible, although it is not gener- 
ically embeddable because it may cause self-intersection 
of the surface. 

To couple the flow to the geometry we must specify how 
the application of a particular Pachner move is triggered 
by the particle content. The rules for fixed geometry in- 
volve particles on a single triangle. The locality of a rule 
which changes the geometry is determined by the local- 
ity of the Pachner moves. The one-to-three, three-to-one 
and two-to-two moves are triggered by the state of one, 
three or two triangles respectively. As in one dimension, 
the constraint of time reversibility is applied in order to 
restrict the set of rules considered [3] . We first recall the 
distinction between invertibility and reversibility. A rule 
is invertible if every state has a unique preimage — given 
the state (particle content plus geometry) one may recon- 
struct the whole unique history leading to that state for 
an invertible rule. For a reversible lattice-gas rule, the 
history of a given state may be generated by an inverse 
rule which can be interpreted as propagation and the 
same collision rules that generate the forward time evo- 
lution. One must recall that the inverse of the product 
of collision and propagation CP, is P~ 1 C~ 1 . 

We choose to apply a one-to-three Pachner move after 
a three particle collision. The particles undergo the three- 
body collision of Figure [3| and then propagate to neigh- 




FIG. 7: New geometry must be created empty so that one 
state does not have two preimages, a problem illustrated in 
this figure. Both figures on the left would produce the figure 
on the right if they were allowed to time evolve. The particles 
on the triangle in the upper left hand corner are therefore re- 
quired to propagate off before the new geometry is created. 
This allows rules with dynamical geometry in which every 
state has a unique preimage, and which are therefore invert- 
ible. 

boring triangles. The Pachner move is applied in the 
subsequent collision. The restriction of time-reversibility 
is satisfied if we create new geometry after the particles 
have propagated from the triangle. That is, we do not 
create new geometry that contains particles. If we cre- 
ated new geometry with particle content, the resulting 
state may have two preimages: one preimage in which 
the geometry is about to be created, and one preimage 
in which particles are about to advect onto existing ge- 
ometry. This problem is illustrated in Figure [7] 




FIG. 8: Creation. A triangle is replaced by a tetrahedron. 
This collision is triggered by a three particle collision, where 
three particles enter a site on every other vector (all odd num- 
bered vectors or all even numbered vectors) such that the 
combined momentum of the three particles is zero. The par- 
ticles then propagate away, and the tetrahedron is formed. 

The rules for addition and deletion are illustrated in 
Figures [8] and [9] When three particles enter a triangle on 
all even numbered vectors or all odd numbered vectors 
and then leave, that triangle is replaced by a tetrahe- 
dron (see Figure [8]) . When three particles propagate off 
a tetrahedron in the same manner, the tetrahedron is 
deleted (see Figure [9]). 

We now determine rules for applying the two-to-two 
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FIG. 9: Deletion. Three triangles forming a tetrahedron are 
replaced by a single triangle. This collision is triggered when 
three particles propagate off the tetrahedron, as shown on 
the left above. This also happens when the particles leave the 
tetrahedron on the empty vectors that point to the surround- 
ing triangles. 



move. This move does not create new geometry and 
so it is straightforward to ensure that the rule is time- 
reversible. This move is triggered by two different states: 
a four-particle state, and a two-particle state, shown in 
Figure [10] 




FIG. 10: The two to two move is triggered by two differ- 
ent states: the two particle state, left, and the four particle 
state, right. The particles remain where they are during these 
moves. 



The rules including dynamical geometry are therefore 
modified from the stochastic FHP rules defined in Fig- 
ures[2] and [3] by the fact that the three body rule of Fig- 
ure [3]is followed by the one-to-three Pachner move shown 
in Figure [8] The stochastic two-body rules remain un- 
changed. The rule set also includes the two-to-two Pach- 
ner moves shown in Figure [10] in which the geometry 
changes by the particle states do not. Naturally, the in- 
clusion of the stochastic two body rule of Figure [2] renders 
the model as a whole irreversible. This could be remedied 
by the addition of a rest particle and the replacement of 
the stochastic two body rule by a deterministic rule as 
described in [BJ. Here we avoid the use of a rest particle 
and retain the stochastic two-body rule. 

The virtue of requiring reversibility of the dynamical 
geometry rules is that the resulting model is fundamen- 
tal, allowing study of the origin of thermodynamics in 
its classical version, and in principal allowing a natural 
generalization to a quantum version. Reversibility also 
constrains the rule space to allow definition of a simple 
and relatively natural model. 
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IV. IMPLEMENTATION 

Implementation of the rules defined in the previous 
section presents several challenges. In this section, we 
describe some of the details of our implementation which 
allow the model to be efficiently simulated without refer- 
ence to embedding space coordinates. 

We first distinguish extrinsic geometry from intrinsic 
geometry. When specifying a triangulation, one can use 
an extrinsic definition of the geometry, or an intrinsic 
definition. An extrinsic definition describes the triangu- 
lation by relating it to an ambient or embedding space. 
For example, a tetrahedron can be defined extrinsically 
by giving the Cartesian coordinates of its vertices in three 
dimensional Euclidean space. 

Geometry can also be defined intrinsically, without ref- 
erence to embedding in some higher dimensional space. 
For example, we can define a tetrahedron intrinsically as 
follows. First, we specify that a triangle is defined by 
three points equidistant from each other. Then, we spec- 
ify that we have four triangles, and that each triangle 
shares exactly one edge with every other triangle. This 



is illustrated in FigurejMj We have defined a tetrahedron 
intrinsically. There was no reference to any coordinate 
system, only reference to parts of the triangulation itself. 

For a lattice gas model defined on an arbitrary triangu- 
lation the flux of particles defines a velocity field. A ve- 
locity vector on a triangulation lives in the tangent space 
to the triangulation at that point. In general, transport 
of tangent vectors on manifolds requires a description 
of the relationship between tangent spaces at different 
points on the manifold. For example, when computing a 
covariant derivative on a Riemannian manifold one must 
consider the variation of coordinate basis vectors with 
position on the manifold. The components of the deriva- 
tives with respect to the coordinates of the basis vectors 
are the Christoffel symbols, which specify the connection 
on the manifold. These quantities are intrinsic: they 
may be computed from the metric without reference to 
any higher dimensional embedding space. 

The implementation of our model contains both in- 
trinsic and extrinsic geometry information. The extrinsic 
(embedding space) information is the set of vertex coor- 
dinates, velocity vectors, and particle coordinates of our 
two dimensional simulation in three dimensional space. 
This is used to produce visualizations of the model, such 
as those shown in Figures [4] and [5} It is possible to imag- 
ine situations in which the extrinsic information is cou- 
pled to the intrinsic model dynamics. For example, a 
membrane embedded in a bulk fluid will have dynamics 
driven in part by the embedding space fluid dynamics. 
We only consider model dynamics defined intrinsically. 
In particular this means that the dynamics remains per- 
fectly well defined even if the triangulation is no longer 
isometrically embedable in three dimensional Euclidean 
space. We do include the possibility that moves which 
would be allowed by the intrinsic dynamics are forbid- 
den conditioned on the embedding, however, for all sim- 



ulations described in this paper these constraints were 
inactive. 

The collision rules are defined locally and must con- 
serve mass, momentum and energy of the particles. We 
wish to apply the same collision rule on every triangle 
expressed in terms of the vector labels. In general, trans- 
lation of a triangle from one location on the triangulation 
to another will induce a transformation of the vector la- 
bels. A reflection of the vector labels through one of the 
symmetry axes of the triangle will change the definition 
of momentum between one triangle and another. Because 
of this we restrict to labelings in which the transforma- 
tion relating the vector labellings of any two triangles is 
one of the three proper rotations of the labeled triangle 
shown in Figure [l] 

The propagation rule moves particles from one triangle 
to another. This operation depends on the transforma- 
tion of the labeling of vectors on going from one triangle 
to another. For each triangle each of the six vectors car- 
ries two pieces of connectivity information which define 
this transformation. Firstly each vector carries a trian- 
gle label which gives the triangle reached by propagation 
along that vector. Secondly each vector carries a vec- 
tor label which determines the vector the particle arrives 
at after propagation. Because the vector labellings of 
any two sites are related by one of three rotations, the 
labelling of two adjacent triangles is determined by the 
image of any one of the vectors. Hence the connectivity 
information redundantly determines the geometry. 



A. Implementing the two-to-two Pachner move 

The two-to-two Pachner move changes the orientation 
of two neighboring triangles but not the number of trian- 
gles. After a collision applying such a move it is necessary 
to change the connectivity information of the surround- 
ing triangles. The move is shown in Figure where 
triangles a and b form a rhombus. During the two-to- 
two move, the four vertices of the rhombus undergo a 
cyclic permutation as the rhombus rotates. The con- 
nectivity between the two triangles involved in the move 
and the surrounding triangles must be updated, and the 
positions of the vertices in embedding space will change 
unless the dihedral angle between a and b is that of the 
icosahedron. The connectivity between triangles a and b 
does not change. 

Triangle pairs which are candidates for the two-to-two 
move are identified before the propagation phase. Trian- 
gles in the appropriate states, for example the state of 



triangle a in Figure 10 are identified. Then, their part- 
ner triangle is examined to see if it in the the state shown 
in triangle b in Figure |10| If it does, the move is per- 
formed; connectivity is redefined with the surrounding 
triangles and the vertices of the triangles are updated. 

As noted above, unless the dihedral angle of the two 
original triangles is that of the tetrahedron, applying the 
two-to-two move will result in a triangulation which is 
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FIG. 11: The two-to-two Pachner move. The pair of triangles 
is rotated in the lattice, so that vertex 1 goes to 2, vertex 2 
goes to 3, vertex 3 goes to 4, and vertex 4 goes to 1. The 
relationship between triangles a and b stays the same; only 
the connectivity between each triangle and the surrounding 
triangles is redefined, along with the vertices of each triangle. 




Lo 



FIG. 12: If the new edge, L, is more than a fraction x different 
from Lo, the two-to-two move is prevented. 



not embeddable in three dimensional Euclidean space. A 
control exists in the code which prevents the change in 
the embedded length of edges from deviating from the 
equilateral value by more than a specified fraction. If Lq 

51 i.e., 



is significantly different from L as shown in Figure 12 



the new edge is significantly different from the new edge, 
we do not allow the two-to-two move to be performed. 
Given a fraction x, < x < 1, we determine whether or 
not the change will be performed via the restriction 



(l-x)Lo < \L\ < (l + x)L 



(1) 



What values of x are relevant? For an initially flat 
triangulation, L — \/3Lq, and for a single tetrahedron 
added to an initially flat triangulation L = \Z2Lq. For 
triangles which meet with the dihedral angle of the icosa- 
hedron L = 4>Lq/2 where <f> is the golden ratio. Hence 
for x < [y2 — 1) with an initially flat triangulation no 
two-to-two moves will be performed. For x < (v3 — 1) 
positive curvature added to an initially flat triangulation 
by a one-to-three move is frozen in place at the new ver- 
tex, as no two-to-two moves may be applied involving any 
face of the tetrahedron. For x < (4>/2 — 1) no two-to-two 
moves may be applied to an icosahedron. 

The effect of these moves is therefore to increase the 
edge lengths of the triangles according the the Euclidean 
metric in three-dimensional embedding space. Either one 
may imagine the triangulation inflated by a scale fac- 
tor and embedded isometrically in a higher dimensional 



space in such a way that the three dimensional embed- 
ding is a projection of this higher dimensional embed- 
ding. In this case the restriction specified by x is to 
triangulations whose projections into three dimensional 
Euclidean space are almost isometric. Alternatively, one 
may regard the three dimensional embedding space itself 
as no longer Euclidean. In this case, x represents a bound 
on the deviation of the metric of the three-dimensional 
embedding space from Euclidean. Note that because x 
specifies a ratio between new and old embedded lengths 
this constraint allows triangle edge lengths grow repeat- 
edly by a succession of geometry changing moves. One 
could also implement a constraint which would bound all 
embedded edge lengths above by an additive constant. 



B. Implementing the one-to-three Pachner move: 
Addition. 



Triangles triggered for addition are marked before 
propagation, and undergo changes in geometry during 
collision. First, the embedding space coordinates of the 
apex of the new tetrahedron are determined. The three 
triangles of the new tetrahedron each have two of the ver- 
tices of the triggered triangle and the third vertex is the 
apex. The move may be regarded as making three copies 
of the original triangle and "rotating" each triangle along 
a different edge so that its free vertex becomes the apex. 
In Figure [13J triangle a has been rotated along the 1 — 2 
edge, triangle b has been rotated along the 0—1 edge, 
and triangle c has been rotated along the — 2 edge. One 
of the three new triangles replaces the original. The con- 
nectivity of the new tetrahedron is set to be that shown 



in Figure 15 where a is the original triangle and b and c 
are the two added triangles. 

The curvature at any triangle vertex is equal to six 
minus the number of triangles meeting at that vertex. If 
the one-to-three Pachner move is implemented without 
restriction, vertices of the triangulation with arbitrarily 
large negative curvature may form. This is because the 
one-to-three move adds a new vertex with positive cur- 
vature and increases the number of triangles meeting at 
each of the original three vertices of the triangle by one. 
In order to allow simulations in which the curvature is 
bounded between ±c we forbid addition of tetrahedra on 
a triangle with any vertex with curvature c. Bounding 
the curvature to be ±1 from the original triangulation is 
equivalent to preventing new tetrahedra from forming on 
existing tetrahedra. 



C. Implementing the 3-1 Pachner move: Deletion. 

The deletion rule depends on the state of three trian- 
gles in a tetrahedral configuration. Tetrahedra are iden- 
tified as sets of triangles whose neighbors are neighbors 
using the intrinsic information - this test uniquely speci- 
fies a tetrahedron (see Figure 14 ) . Once tetrahedra have 
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FIG. 13: The 1 — 3 Pachner move (creation). The triangles 
are labelled x, a, b, and c, and the vertices by 0, 1 and 
2. (m n refers to triangle m's nth vertex). When triangle 
x is triggered for addition, copies of triangle x are rotated 
out of the plane of the paper along each edge. For example, 
to produce triangle b triangle x is copied and rotated along 
the edge joined by the vertices xi and xo, matching vertex 
X2 with the apex. Acting similarly for the other two sides, 
triangle x on the left transitions to the tetrahedron on the 
right. Deletion is the inverse of this process, and intersections 
of common vertices and triangle neighbors are used to identify 
the relevant vertices. The apex is identified as the intersection 
of the vertices of triangles a, b, and c, since the apex is the 
only vertex shared by all three triangles. 



been identified, they are checked to see if their particle 
content makes them a pre-image of deletion. 




FIG. 14: A triangle whose neighbors are neighbors belongs to 
a tetrahedron. This is an intrinsic definition of a tetrahedron. 
On the left there are three triangles sharing one edge each with 
the black triangle. They are the black triangle's neighbors. If 
the blue (top), red (bottom left), and yellow (bottom right) 
triangles in the left figure are also neighbors as designated by 
the arrows above, the three triangles fold into a tetrahedron, 
right. 

Deletion presents a computational issue, in that tri- 
angles must be removed from the list in which they are 
stored. This would change the indices of all triangles, 
requiring a relabeling of the whole triangulation, a com- 
putationally demanding process. It is more efficient to 
disregard the triangles that have been deleted, and place 
them on a dummy list. These obsolete triangles are ig- 
nored whenever the state is updated, and hence relabel- 
ing is avoided. To prevent the list of triangles from ex- 
panding too quickly due to addition and deletion, one of 
the three triangles in a deleted tetrahedron is replaced 
with the single replacement triangle, placing the other 
two on the dummy list. This is the inverse of the one-to- 



three addition move in which the original single triangle 
becomes one of the triangles of the new tetrahedron. 

When a tetrahedron is deleted, three triangles are re- 
placed with one. The triangle with the lowest index is 
retained (in Figure 15 let this be triangle a (right)). This 
triangle replaces the base of the tetrahedron, and so after 
deletion it becomes triangle x in Figure 15 (left). Only 



the coordinates of the apex of this triangle are updated, 
since it rotates about its base (in the case of Figure 15 
the 1 — 2 edge). The vertices involved may be defined 
intrinsically, without reference to their embedding space 
coordinates. The vertex that must be updated is the in- 
tersection of the vertices of all three triangles. The new 
vertex location is the vertex shared by the two triangles 
that are not the replacement triangle and which is not the 



apex. In Figure 15 that is the vertex shared by triangles 



b and c, but not shared by triangle a. The coordinates 
of ao are replaced by the coordinates of bo or vertex Co. 

The algorithm that performs deletion of tetrahedra de- 
pends on the fact that triangles can only be rotated in 
the surface, they can not be flipped. It is convenient to 
define deletion in terms of an involution called inversion. 
When a vector is inverted, the vector is mapped to the 
other vector that occupies the same edge. Referring to 
Figure[TJ the pairs are vectors (0, 1), (2,3) and (4,5). If 
an inversion is performed on vector 3, we get vector 2, 
and so on. This involution is used, together with prop- 
agation along the vectors, to redefine the connectivity 
during deletion using only the intrinsic geometry infor- 
mation 

For example, Figure [15] shows a tetrahedron that will 
undergo deletion and be replaced by triangle a. The con- 
nectivity for vectors as, ao, a.4 and a3 (where a^ for vector 
i of triangle a) must be redefined, as they point to trian- 
gles b and c which will be deleted. After deletion, when 
the tetrahedron is replaced by triangle a, a particle occu- 
pying as will propagate to ko if undisturbed by collision. 
Consider the propagation of a fictitious particle from a 5 
to b.4. After a second propagation this particle would 
end up on triangle c, which is incorrect. Inverting the 
position of the particle, so that it now occupies b 5 and 
allowing the particle to propagate once more takes it to 
ko, which is correct. This was achieved by propagating 
once, inverting, and propagating again. The full set of re- 
labellings given in terms of propagation and inversion are 
shown in Table [Tj Vectors a.4 and are readily identified 
as the vectors attached to the base of the lowest indexed 
triangle of the tetrahedron. They are updated to point to 
J5 and k respectively as shown in Table|TJ The other two 
vectors may be updated by a similar sequence of prop- 
agation and inversion, but it is more straightforward to 
note that ao is updated to ki which is the inversion of 
kg, and a 3 is updated to $4 which is the inversion of J5. 
Hence the update of vectors ao and a 3 is obtained by 
inverting the images of as and a.4, respectively. 



A second useful involution is reflection, which returns 




FIG. 15: When a tetrahedron undergoes deletion, the con- 
nectivity must be redefined. Here, the tetrahedron is being 
replaced by triangle a. 



a<) — > b 5 — -> k — > ki 

p p i 

as c 4 ^> j 5 ^> j4 

p p I 

a4 — ► eg — > C4 — > j 5 

pip 
a 5 — » b4 ^— > b 5 ^— > ko 

pip 



TABLE I: The sequence of relabellings that occur when the 
tetrahedron shown in Figure [l5| is replaced in a three-to-one 
move by the triangle a. P indicates propagation and I indi- 
cates inversion. 



the vector with opposite velocity on the triangle. A re- 
flection on vector 1, for example, returns vector 4. Be- 
cause deletion is allowed when "spectator" particles are 
present on the six edges of the replacement triangle it 
is necessary to update the particle occupancy of the re- 
placement triangle. It is straightforward, in terms of the 
involutions inversion and reflection, to identify the vec- 
tors whose particle occupancy needs to be translated to 
the replacement triangle. 



D. Preventing degenerate triangulations 

In a combinatorial triangulation each triangle is 
uniquely defined by a set of 3 vertexes: it is combinatori- 
ally unique. Triangulations which do not satisfy this cri- 
terion are degenerate. For example, in a degenerate tri- 
angulation two vertices may be connected by more than 
one edge or triangles may share more than one edge. As 
explained below, unrestricted application of the Pachner 
moves can result in degenerate triangulations 

One form of degeneracy occurs when two of a trian- 



gle's neighbors are the same triangle; in other words, two 
triangles share two edges, or only two triangles meet at 
a vertex. Such a feature resembles a "flap" attached to 
the rest of the triangulation. This type of degeneracy is 
avoided by preventing geometry moves whose post-image 
contains a flap. We now consider the effect of the two-to- 
two, one-to-three and three-to-one moves from the point 
of view of avoiding degenerate triangulations. 

Firstly, a flap may be created by the two-to-two Pach- 
ner move. The two to two move increases or decreases 
the number of triangles at a vertex by one. If three tri- 
angles intersect at a point (the apex of a tetrahedron), 
this will become a flap if two of the triangles are replaced 
with one. The application of the two-to-two move to two 
of the triangles of a tetrahedron will therefore result in a 
flap. Preventing two triangles that arc part of the same 
tetrahedron from undergoing a two-to-two move avoids 
this. 

Second, a degenerate triangulation may not be pro- 
duced by the one-to-three move (creation) . Provided the 
initial triangulation is not degenerate any vertex that is 
not at a boundary is shared by at least three triangles. 
The one-to-three move increases the number of triangles 
at each existing vertex by at least one, so this move can 
only create a flap if there were zero triangles to begin 
with. Hence it is not possible for creation (the one-to- 
three move) to result in a degenerate triangulation 

Third, it is possible for the three-to-one move (dele- 
tion) to produce a flap. The three to one move reduces 
the number of triangles at a vertex by one. This move 
can create a flap if the geometry as a whole is a tetra- 
hedron, or if a tetrahedron is attached to a manifold by 
one edge. For example a tetrahedron could be attached 
to the rest of the manifold via a "neck" . To make this 
explicit, take a tetrahedron and label its faces wxyz. Let 
all the faces be connected except for x and y. Now, take 
two adjacent triangles, a and b, in the manifold that are 
not connected to each other, and glue the loose edge of 
y to triangle a, and glue the loose edge of z to triangle b. 
A triangulation with a feature like this is degenerate be- 
cause the two vertices at the join between the tetrahedra 
and the rest of the triangulation are connected by two 
edges. If the tetrahedron that is attached by one edge 
of two triangles were to undergo a deletion of three of 
its faces, a flap would be created. This is illustrated in 
Figure [16] 

This type of structure can be produced from an ini- 
tial geometry which is a tetraspiral - the triangulation 
which results from successive reflections of each vertex of 
a tetrahedron through the opposing face [22] . In this ge- 
ometry, and any geometry composed of tetrahedra shar- 
ing faces, a single three-to-one move results in two tetra- 
hedra connected by a single edge. A second three-to-one 
move will then result in a flap. 

Such degenerate triangulations are prevented from 
forming by two checks. First one prevents degeneracy 
caused by two-to-two moves by checking prior to a two- 
to-two move that both triangles involved do not belong 
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FIG. 16: A tetrahedron attached by an edge to the main 
body of the manifold. If deletion occurs to any triplet of the 
triangles w — x — y — z, a flap will be created. 

to the same tetrahedron. This prevents the formation of 
flaps, and prevents the formation of tetrahedra attached 
by a shared edge. Second, we check that the three-to-one 
move is only applied to tetrahedra which are attached to 
the rest of the triangulation by three edges. It is straight- 
forward to verify that flaps are not produced in our sim- 
ulations by verifying, for example, that the vertex degree 
of all triangles is bounded below by three. 

How does the prevention of such degenerate triangula- 
tions affect the time reversibility of the model? Consider 
the inversion of a sequence of moves involving the cre- 
ation of a degenerate triangulation. Then disallowing 
the formation of degenerate triangulations corresponds 
to disallowing degenerate triangulations in the preimage 
of the inverse rule. If such triangulations are not allowed 
to form, one must naturally forbid them in the initial 
geometry. Disallowing degenerate triangulations in the 
initial condition is sufficient to maintain reversibility. If 
such triangulations are allowed in the initial data, but not 
in the dynamics, invertibility is violated because states 
exist with two preimages: e. g. one where a flap has been 
removed due to addition, and one where it was prevented 
from forming. 



V. GEOMETRY DYNAMICS 



A. Mean field theory 

In this section we consider the average behavior of the 
number of triangles in the model as a function of time. 
First note that deletion is a rare event compared to addi- 
tion. Addition requires exactly three particles in one of 
two configurations of a single triangle. Deletion requires 
exactly three particles in one of two configurations of 
three triangles. If all configurations of a triangle occur 
with equal probability deletion will be less likely than ad- 
dition simply because it requires correlations between the 
states of more than one triangle. Since addition is more 
common than deletion, both by this argument and by ob- 
servation of actual simulations, we construct a mean field 
prediction for the behavior of our system with only addi- 
tion of geometry. Mean field predictions tend to fail fol- 
low dimensional systems. In the one-dimensional case, 
for example, the lattice grew as but a mean field 
model predicted is pQ. It is therefore of interest to de- 
termine the validity of the mean field prediction in two 
dimensions. 

Let TV represent the number of particles and S rep- 
resent the number of triangles. The mean number of 
particles per site is given by 

N 

p = - < p < 6. (2) 

The probability for a site to undergo addition, P + , will 
be proportional to the probability that three sites are 
occupied and three sites are unoccupied. 

P + ^p 3 (l-pf. (3) 

The expected number of triangles which will undergo ad- 
dition is 

(S+) = SP + ^Sp 3 (l-p) 3 . (4) 

The expected change in the number of triangles is given 
by 

AS* = 2(5+) oc Sp 3 {l-p) 3 . (5) 



In this section we study some aspects of the dynamics 
of the geometry degrees of freedom. In these simulations 
the fluid represented by the particles is quiescent - there 
is no forcing applied and because the initial velocities of 
the particles are assigned randomly the average hydro- 
dynamic velocity field will be zero. In one dimension, 
where the only geometrical degree of freedom is the size 
of the lattice, both numerical simulation and calculations 
for particular sets of initial conditions result in an aver- 
age growth of the lattice size of f3 [IH3]. We perform the 
comparable calculations and mean field theory treatment 
of the two dimensional model. In addition, because the 
varying vertex degree of the triangulation represents a 
local geometrical degree of freedom we also present pre- 
liminary results on the distribution of vertex degree on 
the manifold. 



f N 3 iV 4 iV 5 N 6 \ 

AS « (52-- 3 53 + 3 ^ + ^r) ( 6 ) 

In the limit in which creation dominates deletion and 
the number of particles, N, is conserved, the first term 
in the equation above will dominate. Disregarding the 
last three terms, which will become small as the number 
of triangles, S, grows, we convert this to a differential 
equation and solve: 



S oc fs (8) 

The mean field prediction is therefore that the lattice will 
grow asymptotically as is. 
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B. Results 

Four different types of simulations were performed, 
each without restrictions upon the curvature of the mani- 
fold (except those arising from forbidding degenerate tri- 
angulations) or the embedding. All simulations began 
from an initial icosahedral geometry in which each trian- 
gular face was subdivided into 16 triangles by repeatedly 
bisecting the edges. Simulations were performed with 
only creation, only creation and deletion, only creation 
and the two-to-two rule, and a simulation with creation, 
deletion, and the two-to-two rule. While simulations with 
creation but not deletion are not time reversible, they test 
the hypothesis that the geometry dynamics is dominated 
by addition. Thirty realizations were performed for each 
type of simulation with 10 5 timesteps each to determine 
the number of triangles as a function of time. The data 
was fitted to a power law, S(t) = at b , where S(t) rep- 
resents the number of triangles, in the lattice and t is 
the number of timesteps. Fitting the data to the form 
S(t) = a(t — t ) b gave values of to of order one, showing 
that there is only a short transient before the power law 
growth begins, and so fitting the data to S(t) = at b is 
appropriate. 

To evaluate the goodness of the fit x 2 P er degree of 
freedom for the fit function L(t) = at b was computed: 



x 2 (M) 



1 



n — p 



E 



[(S(fy))-m,a,b)] 



(9) 




(a) 




0.335 0.340 
Fit Exponent, b 



where n is the number of data points, p is the number of 
parameters, in this case 2, a and b, and <Ji is the standard 
deviation on the mean for each (Si), where 



(10) 



and j is the number of realizations. 

For each type of simulation an initial fit using Ori- 
gin 7.0 [53] was obtained (using a Levenberg-Marquardt 
method) and an independent error analysis was per- 
formed by computing x 2 hi the a, b plane. The minimum 
value of x 2 found via this method matches that found by 
Origin 7.0. The parameter uncertainties were obtained 
by this x 2 analysis by allowing x 2 to increase by one 
above the minimum. The uncertainties so obtained are 
larger than those given by Origin, presumably because we 
allow a and b to vary independently. In all four types of 
simulation the exponent value is consistent with a power 
law exponent of 1/3, in agreement with the mean field 
prediction. The data and % 2 analysis for the simulations 
with all Pachner moves is shown in Figure[l7j that for the 
simulations with only the three-to-one (addition) Pach- 
ner move is shown in Figure |18[ with all Pachner moves 
except the three-to-one move (deletion) in Figure 19 and 
with all Pachner moves except the two-to-two move in 
Figure [20] 



(b) 



FIG. 17: All Pachner moves allowed. Number of trian- 
gles as a function of time averaged over 30 realizations for 
100000 timesteps each with no restrictions on the curvature 
or embeddability of the triangulation, and all Pachner moves 
utilized. The symbols in a) show the simulation data every 
1000 timesteps and are larger than one standard deviation of 
the mean. The solid line is a fit created in Origin 7.0 with a 
Levenberg-Marquardt method for L(t) = at b to the complete 
data set of 100000 points. Figure b) shows a contour plot 
of x 2 using a sampling of 200 points evenly spaced along the 
range of a and b. The minimum x 2 value lies at x 2 = 0.0319 
at a — 196 and b — 0.33618 in agreement with the fit found by 
Origin 7.0 and the outer-most contour represents a deviation 
of 1.0 from this minimum. The fitted value of the exponent is 
b — 0.33618 ± 0.0065, consistent with a power law exponent 
of 1/3. 



C. Curvature analysis 

Unlike the one-dimensional model it is possible to de- 
fine a curvature variable at each vertex of the triangu- 
lation. As we do not restrict our triangulations in any 
way in the simulations described above it is of interest to 
quantify how curvature is distributed for a typical real- 
ization. We performed four simulations of a single real- 
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(b) 



FIG. 18: Only addition of tetrahedra. Number of trian- 
gles as a function of time averaged over 30 realizations for 
100000 timesteps each. The simulations included only the 
one-to-three Pachner move with no growth or embedding con- 
trol. The symbols in a) show the simulation data every 1000 
timesteps and are larger than one standard deviation of the 
mean. The solid line is a fit created in Origin 7.0 with a 
Levenberg-Marquardt method for L(t) = at b to the complete 
data set of 100000 points. Figure b) shows a contour plot 
of x 2 using a sampling of 200 points evenly spaced along the 
range of a and b. The minimum \ 2 value lies at \ 2 = 1-106 at 
a = 190.864 and b = 0.343 in agreement with the fit found by 
Origin 7.0 and the outer-most contour represents a deviation 
of 1.0 from this minimum. The fitted value of the exponent 
is 0.337 < b — 0.344 < 0.350, consistent with a power law 
exponent of 1/3. 



ization of the type of simulation displayed in Figures |17| 
18|19 20 A histogram of the vertex degree is shown 
in Figure |21| As expected, by allowing unrestricted ad- 



dition of tctrahedran vertices of arbitrarily high degree 
form in the triangulation. However, most of the vertices 
of the triangulation have degree between three and ten. 
While the data shown in Figure [2l] is insufficient to sup- 
port a detailed quantitative analysis of the distribution 
of vertex degree, it appears by inspection to be consis- 




1 ^34 0.336 0.338 0.340 0.342 0.344 0.346 0.348 0.350 
Fit Exponent, b 



FIG. 19: No deletion. Number of triangles as a function of 
time averaged over 30 realizations for 100000 timesteps each. 
The simulations included both two-to-two and one-to-three 
Pachner moves with no growth or embedding control. The 
symbols in a) show the simulation data every 1000 timesteps 
and are larger than one standard deviation of the mean. The 
solid line is a fit of L(t) = at b to the complete data set of 
100000 points. The fit was created in Origin 7.0 using a 
Levenberg-Marquardt method. Figure b) shows a contour 
plot of x 2 using a sampling of 200 points evenly spaced along 
the range of a and b. The minimum x 2 value lies at x 2 = 0.878 
at a = 194.402 and b — 0.343 in agreement with the fit found 
by Origin 7.0. The outer-most contour represents a deviation 
of 1.0 from this minimum. The fitted value of the exponent 
is 0.336 < b = 0.343 < 0.349, consistent with a power law 
exponent of 1/3.. 



tent with an exponential distribution for degrees between 
three and ten. Larger vertex degrees appear to be more 
common than that predicted by this trend below vertex 
degree 10, but there is insufficient data to draw conclu- 
sions here. If we denote the number of tetrahedra added 
to the original geometry N±, and tetrahedra added to 
these N2 and so on, an exponential distribution is con- 
sistent with the ratio of 2Vj to N i+ i being a constant. 

In Figure [22] we display visualizations of the triangu- 
lation for a typical realization. This figure shows a sim- 
ulation in which only the three-to-one Pachner move is 
implemented, resulting in a triangulation which is always 
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FIG. 20: No two-to-two moves allowed. Number of tri- 
angles as a function of time averaged over 30 realizations for 
100000 timesteps each. The simulations included both the 
one-to-three and three-to-one Pachner moves with no two- 
to-two moves and without restrictions on the curvature or 
embeddability of the triangulation. The symbols in a) show 
the simulation data every 1000 timesteps and are larger than 
one standard deviation of the mean. The solid line is a fit of 
L(t) — at b to the complete data set of 100000 points. The 
fit was created in Origin 7.0 using a Levenberg-Marquardt 
method. Figure b) shows a contour plot of y 2 using a sam- 
pling of 200 points evenly spaced along the range of a and 
b. The minimum y 2 value lies at \ 2 = 0.666 at a = 197.060 
and b = 0.335 in agreement with the fit found by Origin 
7.0. The outer-most contour represents a deviation of 1.0 
from this minimum. The fitted value of the exponent is 
0.330 < b — 0.335 < 0.342, consistent with a power law expo- 
nent of 1/3. 



isometrically immersible in three dimensions. This simu- 
lation shows that these triangulations self intersect many 
times. This results in the appearance of many small tri- 
angles in the visualization - in fact these are parts of tri- 
angles which intersect each other. Visualizations of real- 
izations with any combination of rules applied all showed 
this feature. In addition, one observes clusters of added 



FIG. 21: Histogram of number of triangles meeting at each 
vertex. Results of four separate simulations for 100000 
timesteps. The initial geometry was an icosahedron with each 
of its faces subdivided into 16 equilateral triangles. Data is 
shown for one simulation with all Pachner moves implemented 
(Black circles), one simulation with only the three-to-one ad- 
dition move implemented (Red upward pointing triangles), 
one simulation with the one-to-three addition move and the 
two-to-two move but no deletion (Blue stars) and one sim- 
ulation with addition and deletion moves but no two-to-two 
move (Green sideways pointing triangles). 

tctrahedra in all simulations. This may be due to the 
fact that added tetrahedra act as scattering centers for 
particles, and so groupings of tetrahedra will naturally in- 
crease the probability of further geometry-changing col- 
lisions by causing particles to spend longer in a given 
region. 



VI. CONCLUSION 

We have presented the first lattice-gas model with dy- 
namical geometry in two dimensions. Our model is an 
extension of the FHP hydrodynamic two-dimensional lat- 
tice gas model, and the one-dimensional dynamical geom- 
etry lattice gas [IHl]. We have defined and implemented 
rules for dynamical geometry by both Pachner moves. 
For a quiescent fluid on an initially icosahedral geometry 
the number of triangles grows as t 1 ^ 3 for all combinations 
of rules simulated. This is in agreement with a mean field 
prediction, a fact of some interest as mean field predic- 
tions generally fail in low dimensions and in fact fail for 
the one dimensional version of this model [TJ [5] . 

Unlike the one-dimensional case, the flat space limit of 
this model is non-trivial: it is the hydrodynamic FHP lat- 
tice gas. For our model as defined it is therefore possible 
to perform simulations in three regimes. Firstly, the fluid 
may be quiescent, and the geometry dynamical, the limit 
studied in this paper. Secondly, the geometry may be 
fixed and non-trivial, and the fluid driven. This regime 
is relevant for the simulation of surface flows on fixed 
background geometry, such as atmospheric flows on the 
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(a) 




(b) 



FIG. 22: Visualizations of the triangulation for a typical re- 
alization, a) Initial condition for all simulations, b) Typical 
geometry after 100000 time steps for a simulation with only 
the three-to-one Pachner move (addition) implemented, with 
no restriction on the curvature. The faces shown here are in- 
dividual triangles - with only the three-to-one Pachner move 
implemented the embedding is always isometric. The appear- 
ance of many smaller triangles is due to the fact that the 
surface is extensively self-intersected. 



sphere, and in experiments with curved soap films [10] . 
Thirdly, the fluid may be driven and the geometry dy- 
namical, a situation of relevance to surface flows in fluid 
interfaces. Indeed, the equations for surface flow are well 
known, including the case in which the surface is dynam- 
ical mES]. 

Given the model defined here, a natural question to 



pose is: what are the macrodynamical equations of mo- 
tion? In the regime where a non-trivial flow occurs on a 
fixed background geometry, are the relevant fluid equa- 
tions the Navier-Stokes equations on the manifold rep- 
resented by the triangulation? In the regime where the 
geometry is dynamical, does the time evolution of the 
flow coupled to the geometry obey the continuum equa- 
tions of surface flow given in [25] ? 

The principal tool used to obtain the macrodynamical 
equations of a given lattice-gas model is the Chapman 
Enskog expansion |26| . This is an asymptotic expansion 
around a local equilibrium distribution. It is valid in the 
regime that local equilibrium, characterized by a few hy- 
drodynamic fields, is reached rapidly, while global equi- 
libration occurs on longer timescales by hydrodynamic 
processes. Analysis of the model defined in this paper 
requires a new variation of the Chapman Enskog expan- 
sion. 

To treat the model on a fixed, curved surface the Chap- 
man Enskog analysis would need to be extended to arbi- 
trary two-dimensional manifolds. To treat the case where 
the geometry become dynamical it must be possible to 
introduce the geometry degrees of freedom into the Chap- 
man Enskog analysis. One way to do this is to define the 
continuum limit of the triangulation in the same way 
as the continuum limit of the velocity field. That is, 
one considers an average (time, spatial or ensemble) over 
many triangulations of the same surface. The macrody- 
namical equations of surface flow given in [241 125] might 
then arise, for suitably chosen collision rules, from a 
Chapman-Enskog analysis as the slow relaxation of fluid 
plus geometry after a fast relaxation to an equilibrium 
geometry. If such an analysis is valid for the model de- 
fined in this paper, it would also allow simulation of fixed 
geometries via simulation of an ensemble of dynamical 
geometries fluctuating about an average continuum sur- 
face. 

The equilibrium statistical mechanics of two dimen- 
sional triangulated surfaces embedded in three dimen- 
sions has been well studied [ToUTS] . The model defined 
here differs from this body of work in several ways. The 
tethered surfaces studied in |16j have a fixed internal met- 
ric and a Hamiltonian which depends only on extrinsic 
embedding coordinates. The triangulations of our model 
have an intrinsic metric which varies dynamically due 
to application of the Pachner moves. In the terminology 
of [16] this makes our surfaces liquid rather than tethered. 
The object of study of [H] and subsequent work was the 
equilibrium properties of embedded surfaces, here we are 
interested in the non-equilibrium dynamics of surfaces on 
which there is a non-trivial vector field whose dynamics 
is coupled to the intrinsic geometry of the triangulation. 

However, for the case simulated in this paper in which 
the fluid degrees of freedom are quiescent it is interesting 
to compare the typical geometries shown in Figures 22 
with the equilibrium geometries in the crumpled phase of 
random surface models. Two observations are relevant. 
Firstly, as we allow self intersection and do not restrict 
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the embedding our surfaces are phantom surfaces and 
should be compared with random surface models which 
have no extrinsic curvature term in their Hamiltonian. 
The equilibrium geometries of such random surface mod- 
els are crumpled and contain many "spiky" features. 

The typical geometries in our simulations exhibit simi- 
lar features - the high degree vertices shown in Figure [21] 
occur at branching points where many tetrahedra share 
a common vertex. The geometries shown in Figure 22 
also exhibit a concentration of new tetrahedra - showing 
that tetrahedra are added on new tetrahedra more often 
than on the original geometry. This can be explained by 
the fact that the curvature produced by new geometry 
will act as scattering centers for the particles - causing 
particles to spend longer in the vicinity of new geometry, 
where they will then scatter and add further new geme- 
try. This will naturally lead to a branched polymer-like 
structure where tetrahedra are added to tetrahedra and 
particles become trapped on the new branches of the ge- 
ometry. 

Future study of the model should determine whether 
the geometries produced by the model without constraint 
are indeed in the crumpled phase. Simulations in which 
constraints are applied to the local curvature or embed- 
ding of the triangulation may result in geometries closer 
to smooth manifolds and so may be necessary for appli- 
cations in which one aims to simulate a fluid moving on 
a smooth two dimensional surface. 

Fluids, while frequently treated as continua, are in fact 
composed of atoms or molecules. The lattice-gas and 
lattice-Boltzmann methods use the existence of an under- 
lying statistical description of a fluid to realize efficient 
numerical methods for fluid simulation [271 EE] . While 
a discretization of space and time underpins most nu- 
merical methods for field theories, our most fundamental 
current understanding is that spacetime is a continuous 
Lorentzian manifold. 

The idea that, like fluids, spacetime may have under- 



lying discrete substructure occurs repeatedly in specula- 
tive models for quantum gravity. A treatment of classical 
general relativity on polyhedral simplicial complexes was 
first considered by Regge [12]. In the causal dynamical 
triangulations approach the four space-time dimensions 
emerge from an ensemble of simplicial complexes, suit- 
ably constrained by causality [55]. In the causal sets ap- 
proach the underlying Lorentzian manifold is replaced by 
a discrete set of points with a causal (partial) order [3D] • 
In loop quantum gravity geometrical operators such as 
area and length have a discrete spectrum [31 j - In dis- 
crete models of quantum gravity the apparently continu- 
ous classical spacetime emerges at large scales due to the 
smallness of the Planck length. In the more experimen- 
tally accessible world of fluid dynamics, the continuum 
picture is valid because of the largeness of Avagadro's 
number. The discrete model of fluid mechanics on ar- 
bitrary triangulated surfaces presented here provides a 
model in which the question of the emergence of smooth 
manifold-like structures, and associated dynamics of a 
classical field on the manifold, may be studied without 
the numerous conceptual problems of both general rela- 
tivity and quantum mechanics. 
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